# Load packages
library(ggplot2)
library(dplyr)
library(lme4)
library(stringr)
library(sensemakr)
library(stargazer)
library(plm)
# Set seed
set.seed(1996)
#############################
rm(list=setdiff(ls(), c('script', 'scripts', 'log_file')))
#############################
# Load Data
master <- readRDS('incumbent_challengerdf.rds')

# Create IV (WEB) based on challenger position and number of candidates challenging
master$challenger_factor_web <- master$challenger_position_web_1
master$challenger_factor_web[master$candvotes == 1] <- 'No Challenger'
master$challenger_factor_web <- factor(master$challenger_factor_web, c('No Challenger', 'Moderate', 'Extreme'))

# Create IV (CFScore) based on challenger position and number of candidates challenging
master$challenger_factor_cfscore <- master$challenger_position_cfscore_1
master$challenger_factor_cfscore[master$candvotes == 1] <- 'No Challenger' 
master$challenger_factor_cfscore <- factor(master$challenger_factor_cfscore, c('No Challenger', 'Moderate', 'Extreme'))

#############################
# Checking Cross-Pressured Incumbents -- footnote 13
#############################
# See how many races have more than three candidates (2 challengers)
table(master$candnumber >= 3)/nrow(master)

# Creating Cross-Pressure Variable
master$cross_pressure <- NA
# Loop over incumbents
for(i in 1:nrow(master)){
  # If one challenger or fewer, not cross-pressured
  if(master$candnumber[i] <= 2){
    master$cross_pressure[i] <- 0
  }else{ # If more than one challenger, check to see if any qualify as moderate and extreme
    moderate <- sum(str_detect(master[i,str_detect(colnames(master), 'challenger_position_web')], 'Moderate'), na.rm = TRUE)
    extreme <- sum(str_detect(master[i,str_detect(colnames(master), 'challenger_position_web')], 'Extreme'), na.rm = TRUE)
    master$cross_pressure[i] <- ifelse(moderate >= 1 & extreme >= 1, 1, 0)
    
  }
}
# Count number of cross-pressure among multiple challenger races
sum(master$candnumber >= 3 & master$cross_pressure == 1)/sum(master$candnumber >= 3)
# Count number of cross-pressure among all incumbents
table(master$cross_pressure)/nrow(master)

#############################
# Figure 5: Incumbent Challenger Status by Party and Year
#############################
# Create plot df for 2020 and 2018 races with challenger positioning WEB Scores
plot <- master %>%
  filter(year <= 2020) %>%
  filter(!is.na(challenger_factor_web)) %>%
  group_by(year, cand_party, challenger_factor_web) %>%
  summarise(count_chal = n()) %>%
  ungroup()

# Create plot df for 2020 and 2018 races with total challenger counts
plot2 <- master %>%
  filter(year <= 2020) %>%
  filter(!is.na(challenger_factor_web)) %>%
  group_by(year, cand_party) %>%
  summarise(count_total = n()) %>%
  ungroup()
  
plot_df <- merge(plot, plot2)  
plot_df$challenger_prop <- plot_df$count_chal/plot_df$count_total


ggplot(plot_df, aes(x = challenger_factor_web, y = challenger_prop, fill = as.factor(year))) +
  geom_bar(position = "dodge", stat = 'identity') +
  theme_bw() +
  facet_grid(. ~ cand_party) +
  scale_fill_manual(values = c('gray1', 'gray55'), name = 'Election Year') +
  xlab('\nChallenger Status') +
  ylab('Proportion of Incumbents\n') +
  theme(
    # Add these lines to customize the facet_wrap background and text
    strip.background = element_rect(fill = "white"),
    strip.text = element_text(size = 25),
    axis.text = element_text(size = 20),
    axis.title = element_text(size = 25),
    legend.title = element_text(size = 25),
    legend.text = element_text(size = 20),
    legend.position = 'bottom')

ggsave('fg3.tiff', width = 14, height = 8, units = 'in')

# Print counts and percentages for manuscript
plot_df
prop.table(table(master$challenger_factor_web))

#############################
# Table 4: Main Analysis
#############################
# Flip DV for Democrats by multiplying by negative one for WEB Scores
master$model_webscore <- ifelse(master$cand_party == 'Democrat', 
                                master$web_score*-1, 
                                master$web_score)

# Flip DV for Democrats by multiplying by negative one for CFScores
master$model_cfscore <- ifelse(master$cand_party == 'Democrat',
                               master$cfscore*-1, 
                               master$cfscore)

# Coverage in IV
sum(master$candvotes[master$year <= 2020] != 1)
sum(is.na(master$challenger_factor_web[master$year <= 2020]))/sum(master$candvotes[master$year <= 2020] != 1)
sum(is.na(master$challenger_factor_cfscore[master$year <= 2020]))/sum(master$candvotes[master$year <= 2020] != 1)

### Random Effects Models ###
# WEB Score
m1 <- lmer(model_webscore ~ challenger_factor_web + as.factor(year) + (1|FECCandID), 
           data = subset(master, year != 2022))

# CFScore
m2 <- lmer(model_cfscore ~ challenger_factor_cfscore + as.factor(year) + (1|FECCandID), 
           data = subset(master, year != 2022))

### Fixed Effects Models ###
# WEB Score
m3 <- lm(model_webscore ~ challenger_factor_web + as.factor(year) + as.factor(FECCandID), 
           data = subset(master, year != 2022))

# CFScore
m4 <- lm(model_cfscore ~ challenger_factor_cfscore + as.factor(year) +  as.factor(FECCandID), 
           data = subset(master, year != 2022))

# Print stargazer output
stargazer(m3, m1, m4, m2,
          keep = c('challenger_factor_web', 'challenger_factor_cfscore', 'Constant'))



#############################
# Table 5: Sensitivity Analysis
#############################
# Re-orient factor variable with moderate as baseline
master$challenger_factor_web <- factor(master$challenger_factor_web, c('Moderate', 'No Challenger', 'Extreme'))

# Re-estimate model with different challenger baseline 
m3 <- lm(model_webscore ~ challenger_factor_web + as.factor(year) + as.factor(FECCandID), 
         data = subset(master, year != 2022))

# Conduct Sensitivity Analysis with Extreme (ref: Moderate)
sensitivity <- sensemakr(model = m3, 
          treatment = "challenger_factor_webExtreme",
          benchmark_covariates = 'as.factor(FECCandID)H6AZ05083',
          kd = 1:90)

# Print Table output
ovb_minimal_reporting(sensitivity)


#############################
# Table H1: Binary Coding Cross-Pressure
#############################
# Cross-Pressure Binary
master$extreme_web <- NA
# Set extreme to 0 if no challenger or moderate challenger
master$extreme_web[master$challenger_factor_web == 'No Challenger'|
                     master$challenger_factor_web == 'Moderate'] <- 0
# Set extreme to 1 if cross pressured with moderate as the primary challenger
master$extreme_web[master$cross_pressure == 1 & 
                     master$challenger_factor_web == 'Moderate'] <- 1
# Set extreme to 1 if extreme primary challenger
master$extreme_web[master$challenger_factor_web == 'Extreme'] <- 1

# Set moderate to 0 if no challenger or moderate challenger
master$moderate_web <- NA
# Set moderate to 1 if cross pressured with moderate as the primary challenger
master$moderate_web[master$challenger_factor_web == 'No Challenger'|
                     master$challenger_factor_web == 'Extreme'] <- 0
# Set moderate to 1 if cross pressured with extreme as the primary challenger
master$moderate_web[master$cross_pressure == 1 & 
                      master$challenger_factor_web == 'Extreme'] <- 1
# Set moderate to 1 if moderate primary challenger
master$moderate_web[master$challenger_factor_web == 'Moderate'] <- 1

m3 <- lm(model_webscore ~  moderate_web + extreme_web +  as.factor(year) + as.factor(FECCandID), 
         data = subset(master, year != 2022))
stargazer(m3,
          keep = c('extreme_web', 'moderate_web', 'Constant'))

#############################
# Table H2: Flipped Coding Cross Pressure
#############################
master$challenger_factor_web <- factor(master$challenger_factor_web, c('No Challenger', 'Moderate', 'Extreme'))

# Flip-Coding
master$challenger_factor_web_flipped <- NA 
master$challenger_factor_web_flipped[master$cross_pressure == 0 & master$challenger_factor_web == 'Extreme'] <- 'Extreme' 
master$challenger_factor_web_flipped[master$cross_pressure == 0 & master$challenger_factor_web == 'Moderate'] <- 'Moderate' 
master$challenger_factor_web_flipped[master$cross_pressure == 0 & master$challenger_factor_web == 'No Challenger'] <- 'No Challenger' 
master$challenger_factor_web_flipped[master$cross_pressure == 1 & master$challenger_factor_web == 'Extreme'] <- 'Moderate' 
master$challenger_factor_web_flipped[master$cross_pressure == 1 & master$challenger_factor_web == 'Moderate'] <- 'Extreme' 
master$challenger_factor_web_flipped <- factor(master$challenger_factor_web_flipped, c('No Challenger', 'Moderate', 'Extreme'))

table(master$challenger_factor_web_flipped)

m3 <- lm(model_webscore ~ challenger_factor_web_flipped + as.factor(year) + as.factor(FECCandID), 
         data = subset(master, year != 2022))
summary(m3)
stargazer(m3,
          keep = c('challenger_factor_web_flipped', 'Constant'))


#############################
# Table H3: Replicate Results with 2022 Incumbents 
#############################
# Random Effects Models
# WEB Scores
m1 <- lmer(model_webscore ~ challenger_factor_web + as.factor(year) + (1|FECCandID), 
           data = master)

# CFScore
m2 <- lmer(model_cfscore ~ challenger_factor_cfscore + as.factor(year) + (1|FECCandID), 
           data = master)

# Fixed Effects Models
# WEB Score
m3 <- lm(model_webscore ~ challenger_factor_web + as.factor(year) + as.factor(FECCandID), 
         data = master)

# CFScore
m4 <- lm(model_cfscore ~ challenger_factor_cfscore + as.factor(year) +  as.factor(FECCandID), 
         data = master)

stargazer(m3, m1, m4, m2,
          keep = c('challenger_factor_web', 'challenger_factor_cfscore', 'Constant'))

















